!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!! This module contains the implementation of the finite element
!! basis.
!!
!! \n
!!
!! The finite element basis is defined on the canonical simplex, or
!! master element, defined in the corresponding module mod_master_el.
!! It is important to realize that the representation of the basis
!! provided here is inherently connected with the ordering of the
!! nodes and sides defined in mod_master_el. In this module, the
!! master element is denoted by \f$\hat{K}\f$. The type \c t_base
!! defined here is very general and includes, for a chosen degree
!! \f$k\f$,
!! - a vector basis on \f$\hat{K}\f$ (dual basis)
!!   - basis functions \f$\omega_i\f$, \f$i=1,\ldots,m_k\f$
!!   - divergence \f$\nabla\cdot\omega_i\f$
!!   - normal fluxes \f$\omega_i\cdot n_{\partial \hat{K}}\f$
!! - a scalar basis on \f$\hat{K}\f$ (primal basis)
!!   - basis functions \f$\varphi_i\f$, \f$i=1,\ldots,p_k\f$
!!   - gradients \f$\nabla\varphi_i\f$
!! - a scalar basis on \f$\partial\hat{K}\f$ (hybrid basis)
!!   - basis function \f$\eta_i\f$, \f$i=1,\ldots,n_k\f$
!! .
!! Each component of the basis is specified as a symbolic polynomial
!! (see the corresponding module mod_sympoly) as well as through its
!! evaluation at the chosen quadrature points. It is expected that
!! almost all the computations involving the basis can be performed
!! using only the values at the quadrature points; the symbolic
!! representation involves a significant CPU overhead and should be
!! limited to code that is not efficiency critical. The basis
!! description is completed by two quadrature formulas: one for the
!! approximation of \f$\int_{\hat{K}}dv\f$ and one for
!! \f$\int_{\partial \hat{K}}d\sigma\f$. The two quadrature formula
!! are independent from each other and from \f$k\f$, so that it's
!! possible to specify exact quadrature, subintegration or
!! overintegration. It's not necessary to define all the three
!! components of the finite element basis: primal, dual and hybrid.
!! For instance, an implementation of a primal method will only define
!! the primal basis \f$\varphi_i\f$, and the fields corresponding to
!! the dual and hybrid basis will be left unallocated. On the other
!! hand, it's possible to work with more than one t_base object, for
!! instance to build a postprocessing of order \f$k+1\f$ of a solution
!! of order \f$k\f$, or to represent the two variables velocity and
!! pressure in a flow problem. In this cases, of course, one must make
!! sure that the various t_base objects are build with the same
!! quadrature rule.
!!
!! The quadrature rule on \f$\partial\hat{K}\f$ deserves special
!! attention. In fact, as discussed in mod_grid, each side is regarded
!! as the image of the canonical \f$(d-1)\f$-simplex through (at least)
!! three maps: one is the map defined by the ordering of the side
!! vertices and two are the maps defined by the vertex ordering of the
!! neighboring elements. Since these maps are arbitrary, in order to
!! ensure that all of them define the same quadrature rule on the
!! side, the quadrature rule must be invariant with respect to any
!! permutation of the vertices of the \f$(d-1)\f$-simplex or, in other
!! words, it must be consistent with the \f$d!\f$ symmetries of a
!! regular \f$(d-1)\f$-simplex (a simplex with equal edges). For such
!! a quadrature rule, each permutation of the simplex vertices
!! corresponds to a permutation of the integration nodes (and
!! respective weights), and the knowledge of this correspondence
!! between permutations allows us to compute consistently the boundary
!! integrals. The \c t_base type includes a field \c stab to store the
!! permutation table of the boundary quadrature nodes as follows: row
!! \c i of \c stab contains the permutation of the quadrature nodes
!! which corresponds to applying the permutation
!! <tt>me\%pi_tab(i)</tt> to the side vertices.
!<----------------------------------------------------------------------
module mod_base

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   write_octave

 use mod_sympoly, only: &
   mod_sympoly_initialized, &
   t_sympol,       &
   assignment(=),  &
   operator(+),    &
   operator(-),    &
   operator(*),    &
   operator(**),   &
   ev,             &
   pderj,          &
   red_sympoly,    &
   show

 use mod_octave_io_sympoly, only: &
   mod_octave_io_sympoly_initialized, &
   write_octave

 use mod_master_el, only: &
   mod_master_el_initialized, &
   t_me, new_me, clear, &
   write_octave

 use mod_numquad, only: &
   mod_numquad_initialized, &
   get_quad_nodes

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_base_constructor, &
   mod_base_destructor,  &
   mod_base_initialized, &
   t_base,        &
   base, clear,   &
   write_octave

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 !> basis description
 type t_base
   integer :: k !< order (has a different meaning for different basis)
   type(t_me) :: me !< master element \f$\hat{K}\f$
   ! Gaussian quadrature ("nodes" means actually "Gaussian nodes")
   integer :: m  !< # of nodes in \f$\hat{K}\f$
   integer :: ms !< # of nodes for each side \f$s\subset\partial\hat{K}\f$
   integer :: deg  !< exact degree in \f$\hat{K}\f$
   integer :: degs !< exact degree on \f$\partial\hat{K}\f$
   real(wp), allocatable :: xig(:,:)  !< nodes in \f$\hat{K}\f$    (d,m)
   real(wp), allocatable :: wg(:)     !< weights in \f$\hat{K}\f$    (m)
   integer, allocatable :: stab(:,:)  !< side table              (d!,ms)
   real(wp), allocatable :: xigs(:,:) !< nodes in \f$s\f$       (d-1,ms)
   real(wp), allocatable :: wgs(:)    !< weights in \f$s\f$         (ms)
   real(wp), allocatable :: xigb(:,:,:)!< nodes on \f$\partial\hat{K}\f$ (d,ms,d+1)
   real(wp), allocatable :: wgb(:,:)  !< weights on \f$\partial\hat{K}\f$  (ms,d+1)
   ! basis functions
   integer :: mk !< # of basis functions \f$\omega\f$
   type(t_sympol), allocatable :: o_s(:,:) !< \f$\omega\f$ (d,mk)
   type(t_sympol), allocatable :: divo_s(:)!< \f$\nabla\cdot\omega\f$ (mk)
   type(t_sympol), allocatable :: no_s(:,:)!< \f$\omega\cdot n\f$ (mk,d+1)
   real(wp), allocatable :: o(:,:,:) !< \f$\omega(\xi_G)\f$ (d,mk,m)
   real(wp), allocatable :: divo(:,:)!< \f$\nabla\cdot\omega(\xi_G)\f$ (mk,m)
   real(wp), allocatable :: no(:,:,:)!< \f$\omega\cdot n(\xi_{G_b})\f$ (mk,ms,d+1)
   integer :: pk !< # of basis functions \f$\varphi\f$
   type(t_sympol), allocatable :: p_s(:)      !< \f$\varphi\f$ (pk)
   type(t_sympol), allocatable :: gradp_s(:,:)!< \f$\nabla\varphi\f$ (d,pk)
   real(wp), allocatable :: p(:,:)      !< \f$\varphi(\xi_G)\f$ (pk,m)
   real(wp), allocatable :: gradp(:,:,:)!< \f$\nabla\varphi(\xi_G)\f$ (d,pk,m)
   real(wp), allocatable :: pb(:,:,:)   !< \f$\varphi(\xi_{G_b})\f$ (pk,ms,d+1)
   real(wp), allocatable :: gradpb(:,:,:,:)!< \f$\nabla\varphi(\xi_{G_b})\f$ (d,pk,ms,d+1)
   integer :: nk !< # of basis functions \f$\eta\f$
   type(t_sympol), allocatable :: e_s(:) !< \f$\eta\f$ (nk)
   real(wp), allocatable :: e(:,:) !< \f$\eta(\xi_{G_s})\f$ (nk,ms)
 end type t_base

 ! private members

! Module variables

 ! public members
 logical, protected ::               &
   mod_base_initialized = .false.

 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_base'

 interface clear
   module procedure clear_base
 end interface

 interface write_octave
   module procedure write_base_struct
 end interface

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_base_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)  .or. &
       (mod_kinds_initialized.eqv..false.)     .or. &
       (mod_octave_io_initialized.eqv..false.) .or. &
       (mod_sympoly_initialized.eqv..false.)   .or. &
       (mod_numquad_initialized.eqv..false.)   .or. &
       (mod_octave_io_sympoly_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_base_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_base_initialized = .true.
 end subroutine mod_base_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_base_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_base_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_base_initialized = .false.
 end subroutine mod_base_destructor

!-----------------------------------------------------------------------
 
 !> Build a t_base object on the master element \c me.
 !! \details
 !! Except for the master element, all the input variables are
 !! optional. If any of the three finite element basis \c o_s, \c p_s
 !! or \c e_s is omitted, that basis is not generated and the
 !! corresponding fields are left unallocated. If any of the
 !! quadrature formula is omitted, a default is chosen based on the
 !! polynomial order of the basis.
 !!
 !! \note If the arguments \c xig, \c wg and/or \c stab, \c xigs, \c
 !! wgs are present, the arguments \c deg and \c degs are assumed to
 !! be present and to indicate the accuracy of the quadrature
 !! formulas.  In this case, in fact, it is not possible to deduce
 !! such order in this function: it must be known from wherever the
 !! quadrature formula was defined.
 !!
 !! \note The order \c k of the basis is not set in this function,
 !! since it has a different meaning for different combinations of
 !! finite element spaces, such as \f$\mathbb{RT}_k\f$,
 !! \f$\mathbb{DG}_k\f$, \f$\mathbb{BDM}_k\f$ et cetera. Before using
 !! the object returned by \c base, thus, one should set \c k
 !! accordingly.
 !<
 function base(me,o_s,p_s,e_s,deg,degs,xig,wg,stab,xigs,wgs) result(b)
  type(t_me), intent(in) :: me !< the master element \f$\hat{K}\f$
  !> basis \f$\omega\f$, \f$\varphi\f$ and \f$\eta\f$
  type(t_sympol), intent(in), optional :: o_s(:,:), p_s(:), e_s(:)
  !> exactness of the element and side quadrature rules
  integer, intent(in), optional :: deg, degs
  !> quadrature formulas on \f$\hat{K}\f$ and \f$\partial\hat{K}\f$
  integer, intent(in), optional :: stab(:,:)
  real(wp), intent(in), optional :: xig(:,:), wg(:), xigs(:,:), wgs(:)
  type(t_base) :: b

  integer :: &
    ldeg, ldegs, & ! local variables for deg and degs
    is,          & ! side index
    l,           & ! quadrature node index
    i,           & ! basis function index
    h              ! coordinate index
  real(wp) :: bhat(me%d,me%d-1)

  !0) initialization
  b%me = me
  ! used to compute the default values
  ldeg  = 0
  ldegs = 0

  !1) check the vector f.e. space
  if(present(o_s)) then
    b%mk = size(o_s,2)
    allocate(b%o_s(b%me%d,b%mk),b%divo_s(b%mk),b%no_s(b%mk,b%me%d+1))
    b%o_s = red_sympoly( o_s )
    b%divo_s = 0.0_wp
    do h=1,b%me%d
      b%divo_s = red_sympoly( b%divo_s + pderj(b%o_s(h,:),h) )
    enddo
    do is=1,b%me%d+1 ! side loop
      b%no_s(:,is) = 0.0_wp
      do h=1,b%me%d ! dot product
        b%no_s(:,is) = red_sympoly(b%no_s(:,is) + b%me%n(h,is)*b%o_s(h,:))
      enddo
    enddo
    ldeg  = 2*maxval(b%o_s%deg) ! order of the stiffness matrix
    ldegs = ldeg
  else
    b%mk = -1 ! corresponding allocatable fields left unallocated
  endif

  !2) check the scalar f.e. space
  if(present(p_s)) then
    b%pk = size(p_s)
    allocate(b%p_s(b%pk),b%gradp_s(b%me%d,b%pk))
    b%p_s = red_sympoly( p_s )
    do h=1,b%me%d
      b%gradp_s(h,:) = red_sympoly( pderj(b%p_s,h) )
    enddo
    ldeg  = max(ldeg,2*maxval(b%p_s%deg)) ! stiffness matrix order
    ldegs = max(ldeg,ldegs)
  else
    b%pk = -1 ! corresponding allocatable fields left unallocated
  endif
 
  !3) check the hybrid f.e. space
  if(present(e_s)) then
    b%nk = size(e_s)
    allocate(b%e_s(b%nk))
    b%e_s = red_sympoly( e_s )
    if(present(o_s).and.present(p_s)) then
      ldegs = max(maxval(b%o_s%deg),maxval(b%p_s%deg)) + maxval(b%e_s%deg)
    elseif(present(o_s)) then
      ldegs = maxval(b%o_s%deg) + maxval(b%e_s%deg)
    elseif(present(p_s)) then
      ldegs = maxval(b%p_s%deg) + maxval(b%e_s%deg)
    else
      ldegs = 2*maxval(b%e_s%deg)
    endif
  else
    b%nk = -1 ! corresponding allocatable fields left unallocated
  endif

  !4) compute internal quadrature nodes
  if(.not.present(xig)) then ! use deg or ldeg
    if(present(deg)) ldeg = deg
    call get_quad_nodes(b%me%d,ldeg,b%xig,b%wg)
    b%deg = ldeg
  else ! use the input values
    allocate(b%xig(size(xig,1),size(xig,2)))
    b%xig = xig
    allocate(b%wg(size(wg)))
    b%wg = wg
    b%deg = deg
  endif
  b%m = size(b%wg)

  !5) compute boundary quadrature nodes
  if(.not.present(xigs)) then ! use degs or ldegs
    if(present(degs)) ldegs = degs
    call get_quad_nodes(b%me%d-1,ldegs,b%xigs,b%wgs,b%stab)
    b%degs = ldegs
  else ! use the input values
    allocate(b%stab(size(stab,1),size(stab,2)))
    b%stab = stab
    allocate(b%xigs(size(xigs,1),size(xigs,2)))
    b%xigs = xigs
    allocate(b%wgs(size(wgs)))
    b%wgs = wgs
    b%degs = degs
  endif
  b%ms = size(b%wgs)

  !6) complete xigb and wgb
  allocate(b%xigb(b%me%d,b%ms,b%me%d+1),b%wgb(b%ms,b%me%d+1))
  do is=1,b%me%d+1
    ! matrix transformation
    do i=1,b%me%d-1
      bhat(:,i) = b%me%xv(:,b%me%isv(i+1,is))-b%me%xv(:,b%me%isv(1,is))
    enddo
    do l=1,b%ms
      b%xigb(:,l,is) = matmul(bhat,b%xigs(:,l))+b%me%xv(:,b%me%isv(1,is))
    enddo
    b%wgb   (:,is) = b%me%a(is)*b%wgs/sum(b%wgs)
  enddo

  !7) evaluate the basis functions
  if(present(o_s)) then
    allocate( b%o   (b%me%d,b%mk,b%m) , &
              b%divo       (b%mk,b%m) , &
              b%no (b%mk,b%ms,b%me%d+1) )
    do i=1,b%mk ! basis functions
      do h=1,b%me%d
        b%o(h,i,:) = ev(b%o_s(h,i),b%xig)
      enddo
      b%divo(i,:) = ev(b%divo_s(i),b%xig)
      do is=1,b%me%d+1
        b%no(i,:,is) = ev(b%no_s(i,is),b%xigb(:,:,is))
      enddo
    enddo
  endif
  if(present(p_s)) then
    allocate( b%p            (b%pk,b%m) ,         &
              b%gradp (b%me%d,b%pk,b%m) ,         &
              b%pb           (b%pk,b%ms,b%me%d+1),&
              b%gradpb(b%me%d,b%pk,b%ms,b%me%d+1) )
    do i=1,b%pk ! basis functions
      b%p(i,:) = ev(b%p_s(i),b%xig)
      do h=1,b%me%d
        b%gradp(h,i,:) = ev(b%gradp_s(h,i),b%xig)
      enddo
      do is=1,b%me%d+1
        b%pb(i,:,is) = ev(b%p_s(i),b%xigb(:,:,is))
        do h=1,b%me%d
          b%gradpb(h,i,:,is) = ev(b%gradp_s(h,i),b%xigb(:,:,is))
        enddo
      enddo
    enddo
  endif
  if(present(e_s)) then
    allocate(b%e(b%nk,b%ms))
    do i=1,b%nk ! basis functions
      b%e(i,:) = ev(b%e_s(i),b%xigs)
    enddo
  endif

 end function base
 
!-----------------------------------------------------------------------

 pure subroutine clear_base( base )
  type(t_base), intent(out) :: base
   call clear( base%me )
   ! All the other fields are allocatable and are automatically
   ! deallocated
 end subroutine clear_base
 
!-----------------------------------------------------------------------
 
 subroutine write_base_struct(b,var_name,fu)
  integer, intent(in) :: fu
  type(t_base), intent(in) :: b
  character(len=*), intent(in) :: var_name

  integer :: i,j
  real(wp) :: dummy(0,0)
  character(len=*), parameter :: &
    this_sub_name = 'write_base_struct'

   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 30' ! number of fields

   ! field 01 : k
   write(fu,'(a)')      '# name: k'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%k,'<cell-element>',fu)

   ! field 02 : me
   write(fu,'(a)')      '# name: me'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%me,'<cell-element>',fu)

   ! field 03 : m
   write(fu,'(a)')      '# name: m'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%m,'<cell-element>',fu)

   ! field 04 : ms
   write(fu,'(a)')      '# name: ms'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%ms,'<cell-element>',fu)

   ! field 05 : deg
   write(fu,'(a)')      '# name: deg'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%deg,'<cell-element>',fu)

   ! field 06 : degs
   write(fu,'(a)')      '# name: degs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%degs,'<cell-element>',fu)

   ! field 07 : xig
   write(fu,'(a)')      '# name: xig'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%xig)) then
     call write_octave(b%xig,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 08 : wg
   write(fu,'(a)')      '# name: wg'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%wg)) then
     call write_octave(b%wg,'r','<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 09 : stab
   write(fu,'(a)')      '# name: stab'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%stab)) then
     call write_octave(b%stab,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 10 : xigs
   write(fu,'(a)')      '# name: xigs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%xigs)) then
     call write_octave(b%xigs,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 11 : wgs
   write(fu,'(a)')      '# name: wgs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%wgs)) then
     call write_octave(b%wgs,'r','<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 12 : xigb
   write(fu,'(a)')      '# name: xigb'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%xigb)) then
     call write_octave(b%xigb,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 13 : wgb
   write(fu,'(a)')      '# name: wgb'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%wgb)) then
     call write_octave(b%wgb,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 14 : mk
   write(fu,'(a)')      '# name: mk'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%mk,'<cell-element>',fu)

   ! field 15 : o_s
   write(fu,'(a)')      '# name: o_s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   write(fu,'(a)')      '# name: <cell-element>'
   write(fu,'(a)')      '# type: cell'
   if(allocated(b%o_s)) then
     write(fu,'(a,i7)')   '# rows: ',size(b%o_s,1)
     write(fu,'(a,i7)')   '# columns: ',size(b%o_s,2)
     do j=1,size(b%o_s,2)
       do i=1,size(b%o_s,1)
         call write_octave(b%o_s(i,j),'<cell-element>',fu)
       enddo
     enddo
   else
     write(fu,'(a,i7)')   '# rows: ',1
     write(fu,'(a,i7)')   '# columns: ',1
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 16 : divo_s
   write(fu,'(a)')      '# name: divo_s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   write(fu,'(a)')      '# name: <cell-element>'
   write(fu,'(a)')      '# type: cell'
   if(allocated(b%divo_s)) then
     write(fu,'(a,i7)')   '# rows: ',1
     write(fu,'(a,i7)')   '# columns: ',size(b%divo_s)
     do i=1,size(b%divo_s)
       call write_octave(b%divo_s(i),'<cell-element>',fu)
     enddo
   else
     write(fu,'(a,i7)')   '# rows: ',1
     write(fu,'(a,i7)')   '# columns: ',1
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 17 : no_s
   write(fu,'(a)')      '# name: no_s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   write(fu,'(a)')      '# name: <cell-element>'
   write(fu,'(a)')      '# type: cell'
   if(allocated(b%no_s)) then
     write(fu,'(a,i7)')   '# rows: ',size(b%no_s,1)
     write(fu,'(a,i7)')   '# columns: ',size(b%no_s,2)
     do j=1,size(b%no_s,2)
       do i=1,size(b%no_s,1)
         call write_octave(b%no_s(i,j),'<cell-element>',fu)
       enddo
     enddo
   else
     write(fu,'(a,i7)')   '# rows: ',1
     write(fu,'(a,i7)')   '# columns: ',1
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 18 : o
   write(fu,'(a)')      '# name: o'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%o)) then
     call write_octave(b%o,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 19 : divo
   write(fu,'(a)')      '# name: divo'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%divo)) then
     call write_octave(b%divo,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 20 : no
   write(fu,'(a)')      '# name: no'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%no)) then
     call write_octave(b%no,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 21 : pk
   write(fu,'(a)')      '# name: pk'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%pk,'<cell-element>',fu)

   ! field 22 : p_s
   write(fu,'(a)')      '# name: p_s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   write(fu,'(a)')      '# name: <cell-element>'
   write(fu,'(a)')      '# type: cell'
   if(allocated(b%p_s)) then
     write(fu,'(a,i7)')   '# rows: ',1
     write(fu,'(a,i7)')   '# columns: ',size(b%p_s)
     do i=1,size(b%p_s)
       call write_octave(b%p_s(i),'<cell-element>',fu)
     enddo
   else
     write(fu,'(a,i7)')   '# rows: ',1
     write(fu,'(a,i7)')   '# columns: ',1
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 23 : gradp_s
   write(fu,'(a)')      '# name: gradp_s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   write(fu,'(a)')      '# name: <cell-element>'
   write(fu,'(a)')      '# type: cell'
   if(allocated(b%gradp_s)) then
     write(fu,'(a,i7)')   '# rows: ',size(b%gradp_s,1)
     write(fu,'(a,i7)')   '# columns: ',size(b%gradp_s,2)
     do j=1,size(b%gradp_s,2)
       do i=1,size(b%gradp_s,1)
         call write_octave(b%gradp_s(i,j),'<cell-element>',fu)
       enddo
     enddo
   else
     write(fu,'(a,i7)')   '# rows: ',1
     write(fu,'(a,i7)')   '# columns: ',1
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 24 : p
   write(fu,'(a)')      '# name: p'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%p)) then
     call write_octave(b%p,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 25 : gradp
   write(fu,'(a)')      '# name: gradp'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%gradp)) then
     call write_octave(b%gradp,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 26 : pb
   write(fu,'(a)')      '# name: pb'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%pb)) then
     call write_octave(b%pb,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 27 : gradpb
   write(fu,'(a)')      '# name: gradpb'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%gradpb)) then
     call write_octave(b%gradpb,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 28 : nk
   write(fu,'(a)')      '# name: nk'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%nk,'<cell-element>',fu)

   ! field 29 : e_s
   write(fu,'(a)')      '# name: e_s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   write(fu,'(a)')      '# name: <cell-element>'
   write(fu,'(a)')      '# type: cell'
   if(allocated(b%e_s)) then
     write(fu,'(a,i7)')   '# rows: ',1
     write(fu,'(a,i7)')   '# columns: ',size(b%e_s)
     do i=1,size(b%e_s)
       call write_octave(b%e_s(i),'<cell-element>',fu)
     enddo
   else
     write(fu,'(a,i7)')   '# rows: ',1
     write(fu,'(a,i7)')   '# columns: ',1
     call write_octave(dummy,'<cell-element>',fu)
   endif

   ! field 30 : e
   write(fu,'(a)')      '# name: e'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   if(allocated(b%e)) then
     call write_octave(b%e,'<cell-element>',fu)
   else
     call write_octave(dummy,'<cell-element>',fu)
   endif

 end subroutine write_base_struct
 
!-----------------------------------------------------------------------

end module mod_base

